library(tidyverse)
library(survminer)
library(survival)

attach(lung) # attach example data from 'survival' package
Surv(time, status) # create Surv object

fit <- survfit(Surv(time, status) ~ sex, # surv time and status, and grouping factor
  data = lung
)

ggsurvplot(fit,
  data = lung,
  pval = TRUE, # add p value
  risk.table = TRUE, # add risk table
  surv.median.line = "hv", # add median line
  conf.int = TRUE
) # add confidence interval

ggsurvplot(fit,
  data = lung,
  conf.int = TRUE,
  fun = "cumhaz"
) # plot cumulative hazard curve

# this function will add curve of group 'all samples'
ggsurvplot_add_all(fit,
  data = lung,
  palette = "rickandmorty"
) # palette can be changed

ggsurvplot(fit,
  data = lung,
  xlab = "Follow up time (d)", # x axis label
  legend.title = "", # remove legend title
  legend.labs = c("Male", "Female"), # set legend text
  break.x.by = 100, # set tick break on x axis
  pval = TRUE,
  risk.table = TRUE,
  surv.median.line = "hv",
  conf.int = TRUE
)

# with my data of CRC -----------
clinic_meta3 <- read_csv('CRC-I/data/clinic/949CRC-20230918-EngVar-IGHG1+FCGR2B.csv')

survfit(Surv(survival_month, PFS) ~ FCGR2B_I232T, data = clinic_meta3) %>%
  ggsurvplot(data = clinic_meta3,
             xlab = "Follow up time (month)",
             legend.labs = c("II", "IT", "TT"),
             pval = TRUE,
             pval.size = 12,
             break.x.by = 20,
             risk.table = TRUE,
             surv.median.line = "hv",
             title = "CRC patients PFS with FCGR2B-I232T",
             font.legend = 14)

survfit(Surv(survival_month, OS) ~ FCGR2B_I232T, data = clinic_meta3) |>
  ggsurvplot(data = clinic_meta3,
  xlab = "Follow up time (month)",
  pval = TRUE,
  pval.size = 12,
  break.x.by = 20,
  risk.table = TRUE,
  surv.median.line = "hv",
  title = "CRC patients OS with FCGR2b-I232T"
)

# evaluate impact of gross type ----------
type1 <- clinic_meta3 |> filter(gross_type == 1)

survfit(Surv(survival_month, OS) ~ FCGR2B_I232T, data = type1) |>
  ggsurvplot(data = type1,
             xlab = "Follow up time (month)",
             pval = TRUE,
             pval.size = 12,
             break.x.by = 20,
             risk.table = TRUE,
             surv.median.line = "hv",
             title = "CRC patients OS with FCGR2B-I232T"
  )

type2 <- clinic_meta3 |> filter(gross_type == 2)

survfit(Surv(survival_month, OS) ~ FCGR2B_I232T, data = type2) |>
  ggsurvplot(data = type2,
             xlab = "Follow up time (month)",
             pval = TRUE,
             pval.size = 12,
             break.x.by = 20,
             risk.table = TRUE,
             surv.median.line = "hv",
             title = "CRC patients OS with FCGR2B-I232T"
  )

# evaluate combination with IGHG1-G396R -----------
library(data.table)
setDT(crc)

crc$goodgene <- 2
crc[I232T < 3 & G396R == 3, goodgene := 1] # create a new column based on existed values
crc[I232T == 3 & G396R < 3, goodgene := 3]

crcpfsCross <- survfit(Surv(crc$month, crc$PFS) ~ crc$goodgene)

ggsurvplot(crcpfsCross,
  data = crc,
  xlab = "Follow up time (month)",
  legend.labs = c("Survival+/+", "Survival+/-", "Survival-/-"),
  pval = TRUE,
  pval.size = 12,
  break.x.by = 20,
  risk.table = TRUE,
  surv.median.line = "hv",
  title = "CRC patients PFS with FCGR2B-I232T"
)

crcosCross <- survfit(Surv(crc$month, crc$OS) ~ crc$goodgene)

ggsurvplot(crcosCross,
  data = crc,
  xlab = "Follow up time (month)",
  legend.labs = c("Survival++", "Survival+", "Survival-"),
  pval = TRUE,
  pval.size = 12,
  break.x.by = 20,
  risk.table = TRUE,
  surv.median.line = "hv",
  title = "CRC patients OS with FCGR2B-I232T"
)

crc[I232T < 3 & G396R == 3, good2gene := 1]
crc[I232T < 3 & G396R < 3, good2gene := 2] # create a new column based on existed values
crc[I232T == 3 & G396R == 3, good2gene := 3]
crc[I232T == 3 & G396R < 3, good2gene := 4]

crcosCross2 <- survfit(Surv(crc$month, crc$OS) ~ crc$good2gene)

ggsurvplot(crcosCross2,
  data = crc,
  xlab = "Follow up time (month)",
  legend.labs = c("FCGR2B+ IGHG1+",
                  "FCGR2B+ IGHG1-",
                  "FCGR2B- IGHG1+",
                  "FCGR2B- IGHG1-"),
  palette = c('blue', 'green', 'orange', 'red'),
  pval = TRUE,
  pval.size = 12,
  break.x.by = 20,
  risk.table = TRUE,
  surv.median.line = "hv",
  title = "CRC patients OS with FCGR2B-I232T"
)

# combine with MSI variable ----------
clinMerge %>%
  filter((MSI == '0' & SNP1 == 3)|(MSI == '1' & SNP1 == 1) ) ->
  crc

crc_OS_MSI <- survfit(Surv(crc$survival_month, crc$OS.x) ~ crc$MSI + crc$SNP1)

ggsurvplot(crc_OS_MSI,
           data = crc,
           xlab = "Follow up time (month)",
             legend.labs = c("RR + MSI",
                             "GG + MSS"),
           # #palette = c('blue', 'green', 'orange', 'red'),
           pval = TRUE,
           pval.size = 12,
           break.x.by = 20,
           risk.table = 'nrisk_cumcensor',
           title = "CRC patients OS",
           tables.height = 0.35
)

# multivar with gender --------

coxph_model <- coxph(Surv(survival_month, OS) ~ as.numeric(FCGR2B_I232T) + sex, data = clinic_meta3)

## need to set all var here to numeric
ggforest(coxph_model, data = clinic_meta3,fontsize = .8)

g1 <- last_plot()

coxph(Surv(survival_month, PFS) ~ I232T + gender + age, data = pkrm_crc) |>
  broom::tidy(conf.int = TRUE) |>
  mutate(p.value = signif(p.value, 2),
         term = fct_relevel(term, 'I232T', after = Inf)) |>
  plot_forest('PFS')

g2 <- last_plot()

g1 / g2

pkrm_crc |>
  mutate(TT232 = I232T == 3) |>
  coxph(Surv(survival_month, OS) ~ TT232 + gender + age, data = _) |>
  broom::tidy(conf.int = TRUE) |>
  mutate(p.value = signif(p.value, 2),
         term = fct_relevel(term, 'TT232TRUE', after = Inf)) |>
  plot_forest('OS')

g1 <- last_plot()

pkrm_crc |>
  mutate(TT232 = I232T == 3) |>
  coxph(Surv(survival_month, PFS) ~ TT232 + gender + age, data = _) |>
  broom::tidy(conf.int = TRUE) |>
  mutate(p.value = signif(p.value, 2),
         term = fct_relevel(term, 'TT232TRUE', after = Inf)) |>
  plot_forest('PFS')

g2 <- last_plot()

g1 / g2

coxph(Surv(survival_month, OS) ~ G396R + I232T + gender + age, data = pkrm_crc) |>
  broom::tidy(conf.int = TRUE) |>
  mutate(p.value = signif(p.value, 2),
         term = fct_relevel(term, 'G396R', after = Inf)) |>
  plot_forest()

coxph(Surv(survival_month, PFS) ~ G396R + gender + age + I232T, data = pkrm_crc) |>
  broom::tidy(conf.int = TRUE) |>
  mutate(p.value = signif(p.value, 2),
         term = fct_relevel(term, 'G396R', after = Inf)) |>
  plot_forest('PFS')

gender1_crc <- filter(pkrm_crc, gender == 1)
gender0_crc <- filter(pkrm_crc, gender == 0)

survfit(Surv(survival_month, PFS) ~ I232T, data = gender1_crc) |>
  ggsurvplot(pval = TRUE,
             risk.table = 'nrisk_cumcensor',
             title = "Female CRC patients PFS",
             legend.labs = c("II","IT","TT"))

survfit(Surv(survival_month, OS) ~ I232T, data = gender1_crc) |>
  ggsurvplot(pval = TRUE,
             legend.labs = c("II","IT","TT"),
             risk.table = 'nrisk_cumcensor',
             title = "Female CRC patients OS")

survfit(Surv(survival_month, PFS) ~ I232T, data = gender0_crc) |>
  ggsurvplot(pval = TRUE,
             legend.labs = c("II","IT","TT"),
             risk.table = 'nrisk_cumcensor',
             title = "Male CRC patients PFS")

survfit(Surv(survival_month, OS) ~ I232T, data = gender0_crc) |>
  ggsurvplot(pval = TRUE,
             legend.labs = c("II","IT","TT"),
             risk.table = 'nrisk_cumcensor',
             title = "Male CRC patients OS")

# with TCGA dataset -----------
coadSurvival <- read_csv("results/coad+read_survival.csv")

coadOS <- survfit(Surv(crcSurvCluster$days_to_last_follow_up, crcSurvCluster$dead) ~ crcSurvCluster$CD4_CTLA4)

ggsurvplot(coadOS,
  data = crcSurvCluster,
  xlab = "Follow up time (days)",
  # legend.labs = c('II','IT','TT'),
  pval = TRUE,
  pval.size = 12,
  risk.table = TRUE,
  surv.median.line = "hv",
  title = "CRC patients OS with FCGR2b-I232T"
)

ggsurvplot(coadOS,
           data = crcSurvival)
